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Introduction 


Conventional  mammography  is  limited  in  its  sensitivity  for  detecting  subtle  tissue’s 
pathological  changes,  since  the  imaging  relies  on  the  small  differences  in  x-ray  attenuation 
between  the  lesions  and  breast  tissues  of  variable  structure.  As  x-ray  traverses  a  breast,  not 
only  does  its  intensity  get  attenuated,  but  its  phase  also  gets  shifted.  The  amount  of  x-ray 
phase  shift  is  proportional  to  x-ray  wavelength  and  the  ray  integral  of  breast  electron  densities. 
Hence  the  phase  sensitive  x-ray  imaging  has  potential  for  greatly  increasing  x-ray  imaging 
sensitivity  and  specificity  and  reducing  radiation  doses  associated  with  the  imaging. 

The  approaches  for  phase-sensitive  x-ray  imaging  fall  into  three  broad  categories:  the 
crystal  analyzer-based  imaging,  the  phase-gratings  based  interferometric  imaging,  and  the 
inline  phase-sensitive  imaging.  Among  the  three  methods  the  inline  phase-sensitive  imaging 
can  be  implemented  without  any  x-ray  optics  such  as  the  costly  crystal  monochromators  and 
crystal  analyzers  or  the  hard-to-fabricate  expensive  diffractive  gratings,  hence  it  is  especially 
suitable  for  clinical  imaging  implementation.  The  setting  for  inline  phase-sensitive  imaging  is 
very  much  similar  to  conventional  x-ray  imaging,  provided  spatially  coherent  x-ray  illumination 
and  large  subject-detector  distance  are  implemented.  The  phase-contrast  manifests  as  the 
interference  fringes  formed  by  the  phase-shifted  x-rays  undergoing  the  Fresnel  diffraction. 

The  long  term  objective  of  the  project  is  to  develop  a  low-dose  and  quantitative  phase  x- 
ray  imaging  technique  for  facilitating  breast  cancer  detection.  The  Specific  Aims  of  the  project 
of  the  three-year  period  are:  (1)  Develop  prototype  phase-imaging  system  enabling  the  phase 
retrieval,  that  is,  the  reconstruction  of  objects  phase-maps.  The  system  hardware  comprises  a 
micro-focus  tube  operating  at  high  tube  voltages,  a  high  resolution  photostimulable  phosphor 
plate  (CR-plate)  based  detector  system.  The  core  algorithms  for  breast  phase-map 
reconstruction  will  be  developed  to  retrieve  a  breast  phase  map  from  a  single  recorded  image. 
(2)  Validate  the  accuracy  of  the  reconstructed  tissue  projected  electron  densities;  validate  the 
many-fold  radiation  dose  reduction  achieved  with  the  proposed  system;  conduct  subjective 
measurements  to  characterize  the  performance  of  the  proposed  system. 


Body 


In  the  second  year  of  this  project,  as  planned  in  the  Statement  of  Work,  we  performed 
two  tasks:  (A).  Continue  to  develop  the  robust  phase  retrieval  algorithms  for  future  phase 
imaging.  It  is  imperative  to  limit  radiation  doses  involved  in  any  clinical  imaging  applications. 
Hence  it  is  important  to  develop  robust  phase  retrieval  approaches  against  the  noise,  and 
thereby  to  enable  radiation  dose  reduction  in  phase  imaging.  In  order  to  develop  our  new  phase 
retrieval  method  to  its  full  potential,  our  efforts  in  the  Year  2  were  focused  on  refining  the 
attenuation-partition  based  iterative  phase  retrieval  method  and  performing  quantitative  and 
systematic  evaluation  of  its  robustness  and  efficiency  compared  to  the  two  most  widely  used 
phase  retrieval  algorithms  in  the  literature.  (B). Complete  the  system  hardware  integration  for 
future  phase  imaging  with  breast  phantoms. 


(A).  Refinement  and  performance  evaluation  of  the  attenuation-partition  based  iterative 
phase  retrieval  method 
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Introduction:  When  x-ray 
traverses  an  object  x-ray  undergoes 
phase-shift  as  well.  The  differences 
in  x-ray  phase  shifts  between 
different  tissues  are  about  one 
thousand  times  greater  than  the 
difference  in  attenuation.  A  tissue 
phase  map  <p(r)  is  a  map  of  the  ray 

integrals  <p(r)  =  ~(ln/X)  j S(r,s)ds  , 

where  S(f,  s)  denotes  the  tissue 
refractive  index  decrement.  In 
phase  contrast  imaging,  an 
important  task  is  to  retrieve  tissue 
phase  maps.  Phase  retrieval  has 
potential  for  improving  tissue 
contrast  and  achieving  quantitative  tissue  characterization.  Moreover,  phase  retrieval  is 
required  for  phase  sensitive  tomography  and  tomosynthesis  for  removing  artifacts  associated  x- 
ray  refraction  and  diffraction. 

The  mostly  used  phase  retrieval  method  is  based  on  the  Transport  Intensity  of  Equation 
(TIE).  According  to  the  TIE-based  method,  one  acquires  two  images:  one  contact  radiograph, 
and  one  phase  contrast  image.  The  phase  map  <p(r)  is  retrieved  by  using  the  equation 


Fig.1 .  Nylon  air  bubble  Images  acquired  at  40  kVp,  SID  = 
1.75  m.  (left)  Conventional  radiograph,  (right)  Phase- 
contrast  image,  M  =  2.8. 


(p(r)  =  -(2nM  /  AR2 )  V-2 


V ( V-2  (M2l/lin  -  A2q  )) / A20 


[Allen  et  al.  2001],  The  inverse 


Laplacian  operator  V  2  contains  a  singularity  at  the  zero-frequency  in  the  Fourier  domain,  we 
believe  that  it  actually  amplifies  the  noise  in  images  and  may  fail  the  phase  retrieval.  Fig.  1 
shows  a  contact  radiograph  and  a  phase  contrast  image  of  a  nylon  air-bubble  wrap  acquired  at 
40  kVp  from  our  previous  experiment.  We  then  used  the  two  images  in  Fig.  1  for  the  phase 
retrieval.  Fig.  2-(left)  is  the  retrieved  phase  map  with  the  TIE-method  [Yan  et  al.  2008]. 
Apparently  the  retrieval  failed  because  of  relatively  high  noise  levels  in  the  "low  dose" 
projections.  Fig.  2-(middle)  is  the  retrieved  phase  map  with  the  TIE-method  and  the  Tikhonov 
regularization  [Tikhonov  et  al.  1977],  The  Tikhonov  regularization  is  commonly  used 
regularization  method  in  phase  retrievals.  It  essentially  consists  of  replacing  the  inverse 


Fig.  2.  Retrieved  phase  maps  of  the  Nylon  air  bubble  wrap  from  the  acquired  images  in  Fig. 
1.  (left)  Erroneous  phase  map  retrieved  with  the  TIE  method,  (middle)  Erroneous  phase 
map  retrieved  with  the  TIE  method  and  Tikhonov  regularization,  (right)  Phase  map 
retrieved  with  the  prior  knowledge  of  nylon's  attenuation  coefficient  and  refraction  index. 
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Laplacian  operator  V  2  by  vy((V2)2  +a2 ) ,  where  a  is  called  the  Tikhonov  regularization 

parameter.  In  essence  Tikhonov  regularization  seeks  the  minimum-norm,  least  squares  solution 
for  phase  retrieval.  The  retrieval  result  with  Tikhonov  regularization  is  still  unsatisfactory,  as  is 
shown  in  Fig.  (2)-(middle).  In  contrast,  the  image  shown  Fig.  2-(right)  is  the  correctly  retrieved 
phase  map  using  a  different  method  that  utilizes  the  prior  knowledge  of  the  single  material 
(nylon)  attenuation  coefficients  and  its  refraction  indices.  Unfortunately,  the  single  material- 
based  method  does  not  apply  to  general  inhomogeneous  tissues  in  medical  imaging.  These 
findings  show  that  the  TIE-based  phase-retrieval  method  suffers  from  intrinsic  instability  when 
the  noise  presets.  The  robustness  of  phase  retrieval  algorithms  against  noises  in  images  is 
critical  not  only  for  retrieval  accuracies  but  also  for  being  able  to  reduce  radiation  doses 
involved,  since  suppressing  x-ray  quantum  noise  requires  high  radiation  doses  in  imaging.  In 
clinical  imaging  applications  it  is  imperative  to  limit  radiation  doses  involved,  hence  it  is 
important  to  develop  robust  phase  retrieval  approach  for  clinical  applications.  For  general 
inhomogeneous  tissue  samples  encountered  in  medical  imaging  applications,  a  new  phase 
retrieval  method  is  pressingly  needed. 

Method:  In  2008  we  proposed  a  new  phase  retrieval  method:  the  attenuation-partition 
based  iterative  phase  retrieval  method  [Yan  et  al.  2008],  We  have  applied  this  method  to  the 
phase  retrievals  results  for  experimental  images  and  found  satisfactory  results.  We  have 
applied  this  new  algorithm  to  retrieve  the  phase  map  of  a  breast  lumpectomy  specimen  from  its 
contact  radiograph  and  phase  contrast  image,  which  were  acquired  in  our  previous 
experiments.  We  found  that  the  noise  spoiled  the  phase  retrieval  with  the  TIE  method,  and  the 
Tikhonov  regularization  could  not  retrieve  the  phase  map  of  the  specimen  either.  But  the 
attenuation-partition  based  method,  taking  only  10  iterations,  retrieved  the  phase  map  of  the 
specimen  [Yan  et  al.  2008],  However,  only  qualitative  evaluation  of  the  retrieved  phase  map 
was  made  at  that  time 

In  order  to  develop  this  phase  retrieval  method  to  its  full  potential,  in  the  Year  2  we 
refined  our  attenuation-partition  based  iterative  phase  retrieval  method  for  speeding  up  its 

convergence  by  modifying  the  algorithm  coding, 
and  performed  quantitative  and  systematic 
evaluation  of  its  robustness  and  efficiency 
compared  to  the  widely  used  phase  retrieval 
methods  in  the  literature.  In  order  to  simulate  the 
morphological  features  of  breast  tissues,  we  first 
construct  a  digital  breast  specimen  model.  We 
adopted  a  radiograph  of  a  breast  lumpectomy 
specimen  with  pixel  values  rescaled  and  the 
localization  wire  removed  with  a  pixel-value 
interpolation.  The  linear  attenuation  coefficients  and 
electron  densities  for  the  50%  glandular  and  50% 
adipose  breast  tissues  are  computed  from  the 
tissue’s  elemental  composition  and  the  interpolated 
elemental  mass  attenuation  coefficients  from  the 
standard  reference.  Moreover,  each  mass 
attenuation  coefficient  is  further  broken  down  to  two 
components:  one  from  x-ray  photoelectric 
absorption  and  coherent  scattering,  and  another  from  incoherent  scattering.  In  addition,  four 
small  ellipsoids  of  calcium  carbonate  (CaC03)  are  embedded  in  the  specimen  model  to 
simulate  the  breast  microcalcification  in  [Yan  et  al.  2010],  To  simulate  the  noise  encountered  in 
the  practice  of  phase  retrieval,  we  added  into  the  data  the  Poisson-distribution  random  noise, 
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Fig.  3.  Flow  chart  of  the  attenuation- 
partition  based  iterative  algorithm  for 
general  phase  retrieval. 
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and  the  error  caused  by  the  misalignment  between  an  object's  radiograph  and  its  phase- 
contrast  image  acquired  at  a  large  sample-detector  distance.  The  misalignment  could  be 
resulted  from  the  shift  or  tilting  of  the  object  or  detector  between  the  two  projections.  In  the 
algorithm  performance  analysis,  we  systematically  compared  the  performance  of  this  new 
algorithm  with  two  widely  used  phase  retrieval  algorithms  in  the  literature,  namely  the 
Gerchberg-Saxton  (GS)  method  and  the  Transport  of  Intensity  Equation  (TIE)  method  [Allen  et 
al.  2001], 

The  flow  chart  of  our  attenuation-partition  based  iterative  algorithm  is  shown  in  Fig.  3.  In 
this  flow  chart  A2  denotes  the  attenuation-contrast  image,  and  A2KNVc\e  Compton-scattering 
contrast  image,  </>  the  phase  map, /the  phase  contrast  projection  image,  and  the  hypothetical 
phase  contrast  image  formed  by  the  sample  electron  densities  only.  Two  transforms  are 
implemented  in  Fig.  3,  one  is  the  Fresnel  diffraction  transform,  the  other  is  the  phase- 
attenuation  duality  transform  [Yan  et  al.  2010],  Our  idea  underlining  the  algorithm  flow  chart  is: 
use  the  correlation  of  tissue  attenuation  cross-section  with  its  phase  cross-section  to  eliminate 
the  singularity  in  the  phase  retrieval  operator.  This  is  an  extension  to  general  samples  of  the 
phase-attenuation  duality  for  soft  tissues  [Wu  et  al.  2005].  we  first  isolates  the  Compton 
scattering  cross-sections  from  tissue’s  total  attenuation  image,  and  correlate  the  Compton 
scattering  cross-section  with  tissue's  phase  cross-section,  and  retrieve  an  approximate  phase 
map  by  the  phase-attenuation  duality  method  [Wu  et  al.  2005],  which  is  singularity-free  and 
intrinsically  stable.  We  then  iteratively  incorporate  the  rest  attenuation  cross-section  into  phase 


Fig.  4.  Comparison  of  the  TIE  method  and  the  AP-based  method  with  the  noise  added  in  the  simulations, 
(a)  The  true  phase  map.  (b)  Contact  radiograph  with  noise,  (c)  Normalized  phase  contrast  image  with 
noise,  (d)  Phase  map  recovered  with  the  TIE  algorithm  but  without  Tikhonov  regularization,  (e)  Phase  map 
recovered  with  TIE  algorithm  and  with  Tikhonov  regularization;  (f)  Phase  map  recovered  with  the  APBA 
after  10  iteration  steps.  Note  that  the  enhanced  track  in  (c)  is  formed  from  the  phase-sensitive  projection 
of  small  residual  variations  along  the  removed  guide-wire  track  in  the  digital  specimen  model. 
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retrieval  by  repeatedly  correcting  errors,  which  is  calculated  as  the  differences  between  the 
measured  image  intensities  and  the  Fresnel  diffraction  estimates  [Yan  et  al.  2010]. 

Result:  In  the  second  period  of  this  project  we  refined  our  attenuation-partition  based 
iterative  phase  retrieval  method  for  speeding  up  its  convergence.  We  systematically  and 
quantitatively  compared  the  performance  of  this  algorithm  with  other  two  widely  used  phase 
retrieval  algorithms,  namely  the  Gerchberg-Saxton  (GS)  method  and  the  Transport  of  Intensity 
Equation  (TIE)  based  method.  The  systematic  comparison  is  conducted  by  analyzing  phase 
retrieval  performances  with  a  digital  breast  lumpectomy  specimen  model.  We  show  that  the 
proposed  algorithm  converges  faster  than  the  GS  algorithm  in  the  Fresnel  diffraction  regime, 
and  is  much  more  robust  against  image  noise  than  the  TIE  algorithm.  As  an  example,  Fig.  4 
(a)-(f)  show  the  comparison  of  the  TIE  method  and  our  attenuation-partition  based  method  for 
the  digital  breast  lumpectomy  specimen  model.  Our  work  was  published  in  Optics  Express  in 
July  2010  [Yan  et  al.  2010],  The  paper  is  enclosed  in  the  Appendix  of  this  report  for  review. 

(B).  Build  the  prototype  system 


Since  the  project  does  not  support  equipment  purchase,  most  of  the  system  components  are 

those  available  from  our 


subcontractor  Dr.  Liu's 


laboratory  at  University  of 
Oklahoma.  Dr.  Liu  and  his 


group  performed  the  works 
of  integrating  the  hardware 
components  to  a  prototype 
system.  The  system 
hardware  integration  have 
been  continued  in  this 


period.  In  last  period, 
during  the  hardware 
integration,  the  tube 
assembly  malfunction 
occurred,  and  eventually 
the  tube  assembly  and  its 
controller  failed.  After 


tedious  efforts  in  dealing 
with  the  OU  offices  and  the 


vendor  (Hamamatsu 
Photonics)  for  the  repair 
costs  and  repair 

Fig.  5.  A  high  energy  x-ray  phase  contrast  imaging  system  arrangements,  finally  a 

working  source  system 

becomes  available  in  this  period.  The  system  integration  is  completed  in  this  period.  The 
completed  High  Energy  X-Ray  Phase  Contrast  Imaging  System,  as  is  shown  in  Fig.  5,  is  briefly 
described  as  follows. 


The  X-ray  Source  (the  red  box  in  Fig.  5):  The  source  is  a  microfocus  x-ray  source  (Model 
L8121-03,  Hamamatsu  Photonics,  Japan)  was  utilized  for  the  system.  The  x-ray  tube  consists  of 
a  tungsten  target  and  a  Beryllium  output  window.  The  tube  has  a  tungsten-target  and  a  uniform 
circular  focal  spot  of  50[imin  diameter  for  75  W  power  loading,  and  a  focal  spot  of  7  fxmfor  tube 
10W,  and  the  tube  can  operate  at  up  to  150  kVp. 
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The  Image  Detector  (Fig.  5):  Since  tissue  phase-contrast  fringes  recorded  is  proportional  to  the 
Laplacian  and  gradient  of  phase-shifts,  hence  the  detector  should  have  high  sampling 
frequency.  The  image  detection  system  was  a  computed  radiography  system  with 
mammography  plates  (Regius  190,  Konica  Minolta  Medical  Imaging,  Wayne,  New  Jersey)  with 
dimensions  of  24  by  30  centimeters  (cm).  The  mammography  plates  were  designed  with  a  focus 
on  optimizing  the  image  formulation  based  on  low  energies,  which  may  not  provide  optimal 
results  at  high  energies.  The  use  of  general  radiography  plates  may  be  beneficial  instead,  due 
to  the  alternate  design  for  utilization  at  high  energies.  However,  the  CR  system  processes  the 
mammography  plates  with  a  pixel  pitch  of  43.75 pm,  while  general  radiography  plates  are 
processed  with  a  pixel  pitch  of  87.5  |am .  Due  to  the  fact  that  the  mammography  images  are 
produced  with  double  the  spatial  resolution,  they  were  used  for  acquisition  of  the  initial  images. 
However,  the  next  step  in  this  research  is  to  investigate  the  tradeoff  between  the  higher  spatial 
resolution  of  the  mammography  plate  with  the  higher  quantum  efficiency  of  the  general 
radiography  plate,  to  determine  which  of  the  plates  provides  optimal  images  at  high  energies. 

The  Geometric  Configuration  Setting  on  an  Optical  Bench  (Fig.  5):  In  order  to  operate  the 
prototype  in  both  conventional  and  phase  contrast  imaging  modes  for  comparisons,  the  imaging 
and  measurement  components  (the  x-ray  tube,  CR-plate)  are  mounted  on  an  optical  rail,  which 
allows  to  vary  the  source-object  distance  R-i  and  object-detector  distance  R2  as  needed.  It 
allows  a  SID  of  1.83  m.  The  object  to  detector  distance  R2  should  be  sufficiently  large  as  well  to 
allow  phase-shifted  x-rays  to  diffract  as  forming  a  phase-contrast  image.  This  is  because  the 
imaging-geometry  should  optimize  the  phase  contrast  visibility  by  balancing  the  conflict 


Fig.  6.  Image  of  the  Mammography  BR3D  phantom  (140  kVp,  6.4  mAs.) 
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requirements  of  x-ray  spatial  coherence,  large  diffraction  fringes  and  finite  sizes  of  the  focal  spot 
size  and  finite  detector  pitch.  In  a  previous  work  we  found  that  the  phase-contrast  that  could  be 
imaged  is  proportional  to  the  relative  phase-visibility  factor  (RPF)  [Wu  et  al.  2004],  From  the 
estimated  relative  phase-visibility  factors  with  different  geometric  setting  for  this  system,  the 
optimal  magnification  factors  are  found  to  be  in  the  range  of  2  to  2.3. 

Preliminary  Testing:  Phantoms  were  used  for  demonstrating  the  system's  functionality  with  high 
energy  x-rays.  One  of  the  test  images  is  shown  in  Fig.  6.  This  is  an  Image  of  the 
Mammography  BR3D  phantom  (Model  020,  CIRS,  Norfolk,  Virginia),  which  provides  tissue 
equivalent  images  for  qualitative  assessment  of  image  quality.  The  phantom  simulates  a 
realistic  clinical  image,  not  only  through  its  composition,  which  is  100%  adipose  and  glandular 
tissues  blended  together  in  an  approximate  50/50  ratio,  but  also  through  the  5  cm  thickness.  In 
addition,  the  fifth  10  mm  layer  of  the  phantom  is  a  target  slab  imbedded  with  targets  such  as 
specks,  masses  and  fibers.  The  image  as  acquired  with  a  high  energy  phase  contrast  prototype 
system  using  a  Hamamatsu  x-ray  source  with  a  focal  spot  size  of  7  pm  and  a  Konica-Minolta 
Computed  Radiography  System  with  mammography  plates.  The  SID  was  1 .83  m,  and  the  SOD 
was  0.91  m  feet,  which  deliver  a  magnification  factor  of  2  for  the  phase  contrast  images.  The 
image  was  acquired  at  an  experimental  setting  of  140  kV,  6.3  mAs  (71  pA,  90s.) 

In  summary,  the  system  hardware  integration  has  been  completed  in  this  period.  The  tests 
showed  that  the  system  is  functional  for  phase  contrast  image  acquisition  with  high  energy  x- 
rays.  The  efforts  for  technique  optimization  are  under  the  way.  Especially,  we  will  continue  to 
seek  the  ways  to  match  the  CR-plate's  quantum  detection  efficiencies  with  the  high  energy  x- 
rays  as  the  tube  operating  at  120-150  kVp.  One  way  is  to  use  the  general  radiography  CR- 
plates  of  a  pixel  pitch  of  87.5  pm ,  and  investigate  the  tradeoff  between  the  higher  spatial 
resolution  of  the  mammography  plate  with  the  higher  quantum  efficiency  of  the  general 
radiography  plate,  to  determine  which  of  the  plates  provides  optimal  phase  contrast  images  at 
high  energies.  Another  possible  better  solution  is  to  keep  reading  the  general  radiology  CR- 
plates  in  the  high  sampling  mode  of  the  Konica  REGIUS  190  Reader.  But  so  far  the  vendor  is 
resisting  to  give  us  the  access  code  for  the  reader  control  for  using  this  unusual  sampling 
method. 


Key  Research  Accomplishments 


•  It  is  imperative  to  limit  radiation  doses  involved  in  any  clinical  imaging  applications. 
Hence  it  is  important  to  develop  robust  phase  retrieval  approaches  against  the  noise, 
and  thereby  to  enable  radiation  dose  reduction  in  phase  imaging.  The  mostly  used 
phase  retrieval  method  is  based  on  the  Transport  Intensity  of  Equation  (TIE).  We 
showed  that  the  TIE  method  fails  to  provide  robust  phase  retrievals  against  image  noise 
and  the  alignment  errors.  In  this  period  we  refined  the  attenuation-partition  based  phase 
retrieval  method  for  speeding  up  the  iteration  convergence,  and  systematically  and 
quantitatively  compared  its  performance  against  that  of  the  TIE-based  method.  A  digital 
breast  lumpectomy  specimen  model  was  constructed  for  facilitating  the  quantitative 
evaluations.  We  showed  that  the  attenuation-partition  based  phase  retrieval  method 
converges  faster  than  the  GS  method,  and  is  much  more  robust  against  image  noise 
and  the  alignment  errors  than  the  TIE-based  method  does. 
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•  The  system  hardware  integration  has  been  completed  in  this  period.  The  tests  showed 
that  the  system  is  functional  for  phase  contrast  image  acquisition  with  high  energy  x- 
rays.  The  efforts  for  technique  optimization  are  under  the  way. 


Reportable  Outcomes 


In  this  project  period  we  have  published  the  research  results  in  peer-reviewed  journals 
and  presented  the  work  in  an  international  conference,  as  listed  in  the  following. 

Peer-Reviewed  Journal  Article: 

A.  Yan,  X.  Wu,  H.  Liu,  Performance  analysis  of  the  attenuation-partition  based  iterative 
phase  retrieval  algorithm  for  in-line  phase-contrast  imaging,  Optics  Express  18:  16074- 
16089  (2010) 

This  journal  is  the  top-ranking  journal  in  terms  of  the  Impact  Factor  among  all  sixty-four 
journals  in  the  field  of  Optics. 


Published  Abstract  and  Conference  Presentation: 

X.  Wu,  A.  Yan,  H.  Liu,  Improving  robustness  of  phase  retrieval  in  x-ray  phase  contrast 
imaging,  Medical  Physics  37:  3357  (2010) 

X.  Wu,  A.  Yan  ,  H.  Liu,  “Improving  robustness  of  phase  retrieval  in  x-ray  phase  contrast 
imaging,”  Oral  Presentation  at  the  52nd  Annual  Meeting  of  the  American  Association  of 
Physicists  in  Medicine,  July  19,  2010,  Philadelphia,  PA. 


Conclusion 


With  the  support  by  USAMRMC  we  have  successfully  conducted  studies  on  x-ray  phase 
imaging  in  the  second  period  of  this  project.  In  this  period,  we  focused  on  the  development  of  a 
new  phase  retrieval  method  that  is  robust  against  the  noise,  and  thereby  enabling  radiation 
dose  reduction  in  phase  imaging.  We  showed  that  the  mostly  used  phase  retrieval  method  (the 
TIE  method)  fails  to  provide  robust  phase  retrievals  against  image  noise  and  the  alignment 
errors.  In  this  period  we  refined  the  attenuation-partition  based  phase  retrieval  method  for 
speeding  up  the  iteration  convergence,  and  systematically  and  quantitatively  compared  its 
performance  against  that  of  the  widely  used  method.  A  digital  breast  lumpectomy  specimen 
model  was  constructed  for  facilitating  the  quantitative  evaluations.  We  showed  that  the 
proposed  method  is  much  more  robust  against  image  noise  and  the  alignment  errors  than  the 
mostly  used  method  does.  This  new  method  enables  large  radiation  dose  saving  in  phase- 
sensitive  imaging 

In  this  period  we  have  completed  the  system  hardware  integration.  The  tests 
demonstrated  that  the  system  is  functioning  for  phase  contrast  image  acquisition  with  high 
energy  x-rays.  The  efforts  for  technique  optimization  are  under  the  way. 
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In  the  coming  third  period  we  will  continue  the  efforts  in  improving  the  prototype's 
performance.  Especially,  we  will  continue  to  seek  the  ways  to  match  the  CR-plate's  quantum 
detection  efficiencies  with  the  high  energy  x-rays  from  the  tube  operating  at  120-150  kVp.  With 
the  planned  use  of  general  radiography  CR-plates  of  a  pixel  pitch  of  87.5  |am ,  we  will  especially 
need  to  incorporate  the  detector  response  into  the  phase  retrieval  algorithms  for  accurate  phase 
retrieval.  We  will  conduct  phase  imaging  experiments  with  custom-made  breast  phantoms,  and 
will  reconstruct  the  phase  maps  of  the  phantoms.  We  will  experimentally  verify  the  radiation 
dose  saving  with  the  phase  imaging  with  high  energy  x-rays,  and  verify  the  accuracies  of  the 
reconstructed  maps  of  the  projected  electron  densities  of  the  imaged  phantoms. 
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Appendices 
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Abstract:  The  phase  retrieval  is  an  important  task  in  x-ray  phase  contrast 
imaging.  The  robustness  of  phase  retrieval  is  especially  important  for 
potential  medical  imaging  applications  such  as  phase  contrast  mammogra¬ 
phy.  Recently  the  authors  developed  an  iterative  phase  retrieval  algorithm, 
the  attenuation-partition  based  algorithm,  for  the  phase  retrieval  in  inline 
phase-contrast  imaging  [1].  Applied  to  experimental  images,  the  algorithm 
was  proven  to  be  fast  and  robust.  However,  a  quantitative  analysis  of  the  per¬ 
formance  of  this  new  algorithm  is  desirable.  In  this  work,  we  systematically 
compared  the  performance  of  this  algorithm  with  other  two  widely  used 
phase  retrieval  algorithms,  namely  the  Gerchberg-Saxton  (GS)  algorithm 
and  the  Transport  of  Intensity  Equation  (TIE)  algorithm.  The  systematical 
comparison  is  conducted  by  analyzing  phase  retrieval  performances  with 
a  digital  breast  specimen  model.  We  show  that  the  proposed  algorithm 
converges  faster  than  the  GS  algorithm  in  the  Fresnel  diffraction  regime, 
and  is  more  robust  against  image  noise  than  the  TIE  algorithm.  These 
results  suggest  the  significance  of  the  proposed  algorithm  for  future  medical 
applications  with  the  x-ray  phase  contrast  imaging  technique. 

©  2010  Optical  Society  of  America 
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1.  Introduction 

Differing  from  the  conventional  x-ray  imaging,  which  relies  on  the  small  differences  in  x-ray 
attenuation  changes  between  tissues  variable  structure,  inline  phase  contrast  imaging  is  based 
on  the  tissues’  phase-shifts  diffraction  from  the  object  to  the  detector.  Since  x-ray  phase-shift 
differences  between  tissue  and  lesions  are  about  one  thousand  times  larger  than  attenuation 
differences  [2,  3,  4],  x-ray  phase  contrast  imaging  has  the  potential  to  enhance  the  lesion  detec¬ 
tion  sensitivity  and  specificity,  and  reduce  the  radiation  dose  compared  to  conventional  x-ray 
imaging.  In  the  inline  phase  contrast  imaging,  the  diffracted  phase-shifts  form  bright  and  dark 
fringes  at  tissue  boundaries  and  this  bright  and  dark  fringes  are  called  edge  enhancement.  The 
edge  enhancement  relies  on  the  spatial  coherence  of  the  x-ray  source,  the  Laplacian  and  gra¬ 
dients  of  x-ray  phase-shifts  caused  by  the  tissue,  and  the  gradients  of  the  objects  attenuation 
[5,  6,  7,  8,  9,  10,  11,  12,  13,  14,  15],  One  procedure  of  phase  contrast  imaging  is  to  disentangle 
tissue  phase-shifts  from  the  mixed  contrast  mechanism  and  recover  the  phase  maps  from  ac¬ 
quired  phase  contrast  images.  This  procedure  is  called  phase  retrieval.  Phase  retrieval  technique 
plays  a  central  role  in  phase  contrast  x-ray  imaging.  By  means  of  phase  retrieval,  one  can  recon¬ 
struct  a  quantitative  map  of  phase-shifts,  a  phase  image  of  the  imaged  object  [4,  7,  14,  16,  17]. 
The  amount  of  x-ray  phase-shifts  (j)  by  tissues  is  determined  by 

0(?)  =  -  (jf)  re  J Pe(r,z)dz=-  rePe,p(r),  (1) 

where  re  is  the  classical  electron  radius,  h  the  Plank  constant,  c  the  speed  of  light,  E  the  x-ray 
photon  energy,  and  pe  p,  the  integration  of  the  electron  density  pe  over  the  x-ray  path,  is  called 
the  projected  electron  density  [2,  3, 4].  So  a  retrieved  phase  map  is  equivalently  a  map  of  imaged 
object’s  quantitative  projected  electron  densities.  Moreover,  phase  retrieval  is  also  a  necessary 
procedure  for  phase-sensitive  volumetric  imaging,  such  as  the  phase  sensitive  tomography  and 
tomosynthesis,  to  acquire  the  artifact  free  3D  images  of  object  attenuation  coefficients  and 
electron  densities  [8,  15,  16]. 

Phase  retrieval  is  based  on  x-ray  propagation  equation  derived  either  from  the  Fresnel  diffrac¬ 
tion  or  the  Wigner  distribution  based  phase-space  formalism  [5,  18,  9,  19,  20].  To  be  specific, 
let  0(r)  be  the  x-ray  phase-shift  caused  by  the  imaged  object,  and  Ao(r)  the  x-ray  transmis¬ 
sion,  or  the  attenuation-map  of  the  object.  Then  the  Fourier  transform  of  the  x-ray  intensity, 
JF  (/)  ( u )  =  fR 2  7(3 c)  exp  [27 lix-  u]  dx,  at  position  away  from  the  object  with  distance  7?2,  of  the 
monochromatic  point  x-ray  source  starting  at  a  place  away  from  the  object  R\  distance,  can  be 
modeled  by  the  following  [9] 
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where  /m  is  the  incident  x-ray  intensity  at  R\.  X  the  wavelength  of  the  monochromatic  point  x- 
ray  source  and  M  =  (R\  +  Rj) / R\  is  the  geometric  magnification.  When  the  FT-Space  Fresnel 
propagator  kXRju2  /M  <C  1,  Eq.  (2)  can  be  simplified  to  the  Transport  of  Intensity  Equation 
(TIE)  [21,4,  9] 

/(M,  +*)  =  ^  {aJ  (£)  -  ^V.  (l)  J .  (3) 

It  is  worthy  to  note  that  Eq.  (3)  is  valid  only  for  moderate  resolution  images.  For  high  resolution 
images,  i.e.  when  the  FT-Space  Fresnel  propagator  nXRiu2 /M  is  close  or  greater  than  n,  any 
phase  retrieval  algorithms  based  on  Eq.  (3)  need  to  be  substituted  to  Eq.  (2)  [22,  23,  24].  In  this 
paper,  the  algorithms  discussed  are  all  based  on  Eq.  (3),  i.e.  the  moderate  resolution  is  satisfied. 

Numerous  algorithms  have  been  suggested  on  how  to  effectively  recover  the  phase-shift  from 
the  phase  contrast  images.  Among  these,  two  algorithms  are  most  widely  used.  One  is  the  TIE 
algorithm  implemented  by  Allen  and  Oxley  in  [25].  The  other  is  the  GS  algorithm  developed 
first  by  Gerchberg  and  Saxton  in  [26]  and  later  improved  by  Fienup  [27,  28].  These  two  algo¬ 
rithms  have  both  their  advantages  and  disadvantages.  The  TIE  algorithm  is  a  direct  approximate 
method  which  is  fast  but  is  unstable  with  noisy  data;  the  GS  algorithm  on  the  other  hand  is  rel¬ 
atively  more  stable  than  the  TIE  algorithm[25]  but  the  converging  speed  is  slow  especially  for 
the  field  of  medical  imaging.  In  [1],  the  authors  developed  an  Attenuation-Partition  Based  Al¬ 
gorithm  (APBA)  based  on  the  phase-attenuation  duality  property  of  soft  tissues  under  higher 
x-ray  energy.  This  algorithm  is  fast  and  stable  for  potential  medical  imaging.  We  compared  the 
performance  of  this  algorithm  with  the  TIE  algorithm  for  two  groups  of  data  under  the  condi¬ 
tion  of  medical  imaging  in  [1],  including  the  phase  retrieval  from  phase-contrast  images  of  a 
breast  lumpectomy  specimen.  In  this  paper,  we  will  make  a  systematic  analysis  about  this  al¬ 
gorithm  and  compare  its  performance  with  the  one  of  the  GS  algorithm  and  the  TIE  algorithm 
with  simulated  data. 

The  paper  is  organized  as  follows.  In  Section  2,  we  first  summarize  the  attenuation-partition 
based  algorithm  (APBA),  which  is  motivated  by  our  observation  of  the  phase-attenuation 
duality [14].  Then  we  give  a  measure,  called  total  variation,  used  to  evaluate  the  closeness  of 
two  image  data.  This  measure  is  used  as  a  quantitative  standard  in  comparing  the  performance 
between  different  algorithms  in  the  following  section.  In  Section  3,  we  first  develop  a  breast 
specimen  model  which  can  reflect  the  attenuation  and  phase  changes  with  respect  to  the  x-ray 
energy  change  (Section3.1),  and  then  compare  the  performance  of  the  algorithm  with  the  GS 
algorithm  (Section  3.2)  and  the  TIE  algorithm  (Section  3.3).  Finally,  we  conclude  this  paper  in 
Section  4. 


2.  The  attenuation-partition  based  algorithm  (APBA)  and  an  image  accuracy  measure 

2.1.  The  Attenuation-Partition  Based  Algorithm 

The  attenuation-partition  based  algorithm  (APBA)  is  a  recently  developed  iterative  algo¬ 
rithm  for  phase  retrieval [  1],  It  was  derived  from  our  previous  notion  of  the  phase-attenuation 
duality[14],  and  it  takes  advantage  of  the  correlation  between  the  attenuation  and  phase-shift  for 
phase  retrieval.  As  is  well  known,  tissue’s  attenuation  change  Ao  in  the  diagnostic  x-ray  imag¬ 
ing  arose  from  three  x-ray  and  tissue  interactions:  the  photoelectric  absorption,  the  coherent 
scattering,  and  the  incoherent  scattering [29,  30,  9,  14].  However,  among  the  three  interactions, 
the  attenuation  caused  by  incoherent  scattering  Akn,  which  is  dominated  by  the  x-ray  Compton 
scattering,  deserves  a  special  attention.  This  is  because  both  of  Akn  and  the  x-ray  phase-shift  <j> 
are  determined  by  the  tissue’s  projected  electron  density; 


Akn(?)  =  exp 


0(c)  —  Arepep, 


(4) 
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where  A  is  the  x-ray  wavelength,  re  the  classical  electron  radius,  pe  p  the  projected  electron 
density  as  defined  in  Eq.  (1),  and  o^n  is  the  total  cross  section  for  x-ray  Compton  scattering 
with  a  free  electron: 


t^KN  (Ephoton  )  —  A?Tre 


1+1  [2(i+i) 

tj2  [  l+2rj  77 


+  —  log(l  +277) 


1  +  377  1 
(I+277)2/’ 


(5) 


with  rj  =  Ephoton/mec2.  Here  we  denote  the  photon  energy  of  the  primary  x-ray  beam  by  Ephoton 
and  mec2  is  the  rest  electron  energy.  Eq.(4)  suggests  clearly  that  the  x-ray  attenuation  and  phase- 
shift  by  tissue  has  certain  correlation.  Our  idea  is  to  utilize  this  correlation  for  facilitating  the 
phase  retrieval. 

Of  course,  the  extent  of  this  correlation  between  phase  and  attenuation  depends  on  the  x- 
rays  photon  energy  as  well  as  the  tissues  physical  composition.  For  example,  for  light  elements 
with  atomic  numbers  Z  <  10,  x-ray  attenuation  is  dominated  by  the  Compton  scattering  for 
x-rays  of  60  keV  or  higher,  i.e.  +0  ~  ,4  kn  [  1 4] .  We  call  this  relationship  between  phase-shift 
and  attenuation  the  phase-attenuation  duality.  The  phase-attenuation  duality  can  be  used  for 
phase  retrieval  as  follows.  Consider  a  phase  contrast  imaging  setting  with  a  point  source  of 
wavelength  A.  The  object  is  at  a  distance  R\  from  the  source.  We  denote  AS  as  the  distance  from 
object  to  detector  plane,  M  =  (R\  +  R2) /R 1  the  geometric  magnification  factor,  Im  the  entrance 
x-ray  intensity  at  R\ ,  and  Id(jd)  the  x-ray  intensity  at  the  detector  plane.  For  convenience, 
we  denote  I  =  M 2  ■  Iq (jD)/kn  as  the  normalized  intensity  of  phase-contrast  image.  When  the 
phase-attenuation  duality  holds,  the  phase  map  (j)  (r)  can  be  robustly  retrieved  from  just  a  single 
projection  image  [14]: 


A|N(r)  =  S)(7)  =  #-1 


(  ^00  \ 

y  1  +  47T2Au2  J  ’ 


(6) 


where 

l=^Rb-—'  (7) 

2nM  Okn 

and  A),  for  sake  of  convenience,  is  called  the  “duality  transform”  acting  on  the  normalized 
image  I. 

In  general  imaging  cases,  such  as  with  low  energy  x-rays  or  an  object  contains  calcified 
tissues  such  as  calcification,  this  phase-attenuation  duality  does  not  hold.  However,  we  can  still 
factor  out  tissue’s  total  attenuation  Ao  as 


A0(?)  =  Akn(?)  •Ape)Coh(r),  (8) 

where  we  denote  the  attenuation  caused  by  photoelectric  absorption  and  coherent  scattering 
by  APe,coh-  Strictly  speaking,  Okn  is  only  Compton  scattering  cross-section,  it  may  be  slightly 
different  from  the  incoherent  scattering  cross-section  for  high-Z  elements.  This  is  because  while 
Compton  scattering  is  x-ray  scattering  from  the  free  electrons,  the  incoherent  scattering  is  that 
from  the  bound  atomic  electrons [31].  So  when  we  factor  Ao  =  Akn  'Ape.coh  (Eq.  (8))  we  actually 
factor  the  small  residual  binding  effect  of  atomic  electrons  into  ApeiCOh-  With  this  understanding, 
Eqs.(4)  and  (8)  are  rigorously  valid.  The  notion  of  Eq.  (6)  and  Eq.  (8)  led  us  to  the  development 
of  the  attenuation-partition  based  algorithm  [1],  While  the  derivations  and  the  algorithm  details 
of  this  method  can  be  found  from  Ref.  [1],  a  brief  description  of  the  method  is  as  follows. 
Our  goal  is  to  retrieve  the  phase  map  from  the  two  normalized  images:  one  is  the  object’s 
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attenuation  image  A5  acquired  with  Ri  =  0,  and  the  other  is  the  acquired  phase  contrast  image 
I  =  M 2  ■  //j(r'/j)  //in  with  Ri  >  0.  Employing  the  acquired  phase  contrast  image  /  and  the  duality 
transform  Eq.  (6),  we  will  first  obtain  an  estimate  for  the  attenuation-component  Akn  and  phase 
map  (j).  We  then  rewrite  Eq.  (8)  as 

^0=^KN  —  8a,  8A  =  Akn(1  —  Ape.coh),  (9) 

and  find  the  correction  term  (5,4  using  the  estimate  of  Akn-  We  then  employing  the  Fresnel 
Diffraction  transform  (as  defined  in  Eq.  (11))  to  transport  the  wavefield  8Ae'^  from  /\’|  to  AS. 
We  can  find  81  =  (8Ael&)  |2,  which  is  the  difference  between  phase  contrast  image  /  and 

the  “duality-only”  counterpart  /kn  =  |3r  (Akn<?!<^)|~-  With  the  corrected  “duality-only”  im¬ 
age  intensity  /kn  =  [\Tl  +  \/8 /  j  we  can  start  a  new  round  of  iterations  by  repeating  above 
procedure.  For  a  rigorous  analysis  of  the  iterative  algorithm  and  its  convergence  interesting 
readers  are  referred  to  [1].  Note  that  the  equation  \/l  =  \JIkn  —  \/<5/  is  generally  valid,  since 
it  is  actually  a  result  of  x-ray  Fresnel  diffraction  and  extremely  smallness  of  hard  x-ray  wave¬ 
length  compared  to  finest  resolution  achievable  in  the  phase  contrast  imaging.  While  interesting 
readers  can  find  the  mathematical  proof  of  this  equation  in  [1],  an  intuitive  explanation  of  this 
formula  comes  from  the  x-ray  propagation.  In  such  a  wave  propagation  process,  or  the  so-called 
semiclassical  wave  propagation,  the  phase  of  a  wave  field  evolves  essentially  according  to  the 
free-space  Hamilton- Jacobi  equation  from  its  initial  phase  value.  So  if  we  denote  the  Fresnel 
free  propagation  as  a  Fresnel  transform  gr  acting  on  the  initial  wavefield,  therefore  the  two 
resulted  wavefields  (Akn  exp [;</)] )  and  (5Aexp[i0])  would  have  the  same  resultant  phases, 
so  the  resultant  amplitude  from  the  two-wave  superposition  is  given  as  \Jl  =  \/Ikn  —  Vo7. 

The  above  iterative  procedure  can  be  summarized  in  flow  chart  Fig.  1  and  the  Algorithm 
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Fig.  1.  Flow  chart  of  APBA 


Algorithm.  In  an  imaging  experiment,  assume  we  have  performed  two  imaging  measurements. 
One  image  (radiograph)  is  the  attenuation  image  Aq  acquired  at  SID=  R\,  another  is  a  nor¬ 
malized  phase-contrast  image  I  =  M 2  ■  Id(/d) /ha  acquired  at  SID=  R \  +  Rz-  With  Aq  and  I  as 

well  as  the  initial  81,  usually  0,  as  known  inputs,  we  first  assume  /kn  =  (f/l  +  V8l^j  ■  Then 
(1).  Compute  Akn  =  \/5)  (/kn)  and  (j)  from  the  duality  transform  Eq.  (6). 
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(2).  Compute  S A  from 


SA  =  Akn  •  (1  —  P),  P  =  Aq/Akm. 


(10) 


Equations  (9)  and  (10)  are  in  fact  the  same  equations.  The  advantage  of  Eq.  (10)  over 
(9)  is  that  we  can  set  a  threshold  for  P.  We  know  P  =  Ape,coh  in  the  ideal  case  and 
Ape.coh  is  bounded  between  1  and  Aq.  The  computation  rounding  error  or  the  presence  of 
measured  noise  in  the  acquired  data  can  make  the  value  of  P  over  pass  these  bounds  in 
the  iterative  computations.  By  setting  a  threshold  upper  bound  ubd  =  1  and  lower  bound 
lbd  =  min(Ao),  the  minimum  value  ofAo,  to  P  in  the  iterative  computations,  we  can  make 
the  algorithm  more  stable.  Moreover,  if  we  know  a  better  lbd  for  Ape,coh,  other  than  the 
minimum  ofAo,  the  converging  speed  of  the  algorithm  can  be  greatly  improved. 

(3).  Compute  81  by  Fresnel  propagate  SAe11^  from  position  R\  to  R2:  Si  =  \fjx(8Ae'^\f  with 
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(4) .  Compute  /kn 


Go  to  (1 )  for  next  iteration. 


The  number  of  iterations  or  an  accuracy  measure  can  be  used  to  determine  when  to  exit  the 
program:  assuming  | -|  is  some  kind  of  norm,  that  can  effectively  reflect  the  accuracy  of  the 
retrieved  data  as  an  image,  if  ||54+i  —  5/*||  is  less  than  a  predefined  threshold  value  /j,  or  the 
iteration  step  exceeds  a  predefined  maximum  number  of  iteration  steps,  then  0  is  the  retrieved 
phase  and  the  iteration  stops;  otherwise  further  iteration  is  needed. 

An  appropriate  image  accuracy  measure  should  be  a  measure  that  can  effectively  reflect 
the  accuracy  of  the  data  AND  at  the  same  time  correctly  reflect  the  visual  perception  of  the 
data  as  an  image,  since  an  image’s  visual  perception  is  crucial  for  diagnostic  radiology.  In  the 
next  subsection  the  authors  suggest  a  measure  which  can  be  employed  as  an  appropriate  image 
accuracy  measure,  which  was  first  investigated  by  Rudin  in  [32], 


2.2.  An  Image  Accuracy  Measure 


(a)  ( b )  (c) 

Fig.  2.  Example  “Lena”  images  used  in  measuring  the  closeness  between  images,  (a)  <j)\ , 
(b)  02,  (c)  03 


A  continuous  signal  is  generally  represented  as  a  function  of  vector  variables:  /( r).  A  sam¬ 
pled  signal  will  be  represented  as  a  one-  (or  higher)  dimensional  sequence  of  real  numbers.  In 
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this  paper,  we  will  denote  the  continuous  two  dimensional  image  as  an  intensity  function  of 
two  dimensional  variables,  such  as  or  /( ?),  r  =  (x,y).  The  sampled  2-D  image  will  be 

represented  by  i  =  1 , 2,  •  •  -  ,m,  j  =  1 , 2,  -  •  •  ,n.  In  practice,  to  estimate  a  true  signal  in 

noise,  the  most  frequently  used  methods  are  based  on  the  least  squares  criteria  and  thus  an  intu¬ 
itive  measure  for  the  closeness  of  two  image  functions  /  and  g  is,  similar  to  statistical  standard 
deviation,  based  on: 


std (g,f)  :=  1 


fn(g(f)-f(f)-n)2df 


1/2 


V(0) 


(12) 


where  £2  is  the  finite  domain  of  image  functions  /  and  g,  V  (£2)  represents  the  area  of  domain 
£2,  and  fi  =  fn(g  —  f)dr/V(£2)  is  the  statistical  mean  value  of  g  —  f.  In  [32],  L.  I.  Rudin 
investigated  the  relation  of  edge  formation  of  the  2-D  digital  image  and  the  singularities  of 
the  image  function  and  pointed  out  that  the  image  intensity  function  belongs  to  the  space  of 
functions  of  bounded  total  variation.  Rudin  et  al.[33]  pointed  out  that  the  proper  norm  for 
images  is  the  total  variation  (TV)  norm,  which  is  the  L\  norm  of  the  gradient  of  the  image 
function,  and  not  the  Li  norm.  For  two  image  functions  /  and  g,  we  define  the  TV  norm  of 
g  —  f  as  the  closeness  of  the  two  images: 


T  V(g,f) 


V(£2) 


dy 


V(£2) 


(13) 


where  V  is  the  gradient  operator.  Since  std  is  the  form  of  Li  norm,  it  is  not  suitable  as  a 
measure  to  represent  the  closeness  between  two  images.  For  example,  for  the  three  “Lena” 
image  functions  shown  in  Fig.  2,  the  std  measures  of  (f)\  —  02  and  01  —  03  are  std(0i ,  0?)  =  0.206, 
std(0i ,  03)  =  0.472  respectively.  But  the  image  0 3  is  much  more  closer  to  0 1  than  0 2  does  in 
visual  perception.  This  is  because  the  digital  representation  of  an  image  depends  not  only  on 
the  pixel  values,  it  also  depends  on  the  contrast  changes  between  neighbor  pixels.  An  image’s 
visual  perception  is  crucial  for  diagnostic  radiology.  This  contrast  changes  between  neighbor 
pixels  can  better  be  represented  by  the  gradient  changes  of  the  image  functions.  For  example, 
the  TV  measures  of  0i  —  0 2  and  0 1  —  0 3  are  TV(0i,02)  =  0.0247  and  TV(0i,03)  =  0.0138 
respectively,  more  appropriate  in  reflecting  the  visual  perception.  In  this  paper,  we  will  use  the 
TV  norm  (13)  as  the  measure  for  closeness  between  two  compared  image  functions. 

Note  that  the  TV  norm  between  two  image  functions,  say  g  and  /,  equals  0  if  and  only 
if  g  differs  from  /  by  a  constant.  This  feature  does  not  affect  its  appropriateness  for  phase 
retrieval  since  it  is  well  known  that  the  recovered  phase  0  is  unique  up  to  a  constant  with  given 
information  about  the  attenuation-map  ,4q  and  the  phase-contrast  intensity-map  I. 


3.  Simulation  Tests 

In  order  to  investigate  the  performance  of  the  algorithm  constructed  above,  we  perform  com¬ 
puter  simulations  in  this  section.  In  Section  3.1,  we  first  construct  a  breast  tissue  model  that 
represents  the  phase-shifts  and  attenuation  of  breast  tissues  and  embedded  microcalcifications 
for  different  x-ray  energies.  In  our  simulation  tests,  the  distances  of  source  point  to  object,  Ri, 
and  object  to  detector,  R2,  are  set  to  26  inches  (0.66  m)  respectively.  In  this  way  the  magni¬ 
fication  factor  M  =  (R\  +  Ry) / R\  is  equal  to  2.  For  convenience,  the  incident  x-ray  intensity 
/in  at  R\  is  set  to  M 2  (4).  We  compare  the  performance  of  our  algorithm,  APBA,  with  two 
other  widely  used  algorithms:  the  Gerchberg-Saxton  (GS)  algorithm  (Section  3.2)  and  the  TIE 
algorithm  (Section  3.3). 
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3.1.  A  breast  specimen  model 

The  tissue’s  phase-shift  and  attenuation  are  determined  not  only  by  the  tissue’s  physical  com¬ 
position,  but  also  by  the  x-ray  energy.  With  different  x-ray  energy,  the  same  tissue  has  different 
phase-shift  and  attenuation  change.  Simulation  models  used  in  literature  often  do  not  incorpo¬ 
rate  these  changes.  In  this  subsection,  we  construct  a  breast  specimen  model  that  can  represent 
tissue  attenuation  and  phase  shifts  according  to  employed  x-ray  energies  as  well  as  tissue’s 
compositions. 


(a)  (b)  (c) 

Fig.  3.  Image  manifest  of  (a)  A^e  coh,  (b)  Aj^  and  (c)  Ag  when  x-ray  energy  equals  35.5  keV. 


Fig.  4.  Profiles  of  A^e  coh,  the  solid  lines,  and  Aj^,  the  dotted  lines,  when  x-ray  energy 
equals  (a)  18.5  keV,  (b)  35.5  keV  and  (c)  59.5  keV. 


In  our  model  the  tissue  has  two  physical  compositions:  the  50%  glandular-50%  adipose 
breast  tissue  and  the  microcalcifications.  In  order  to  simulate  the  morphological  aspects  of 
breast  tissues,  we  adopted  a  radiograph  of  a  breast  lumpectomy  specimen  with  pixel  values 
rescaled  and  the  metal  localization  wire  removed  by  replacing  the  pixel  value  at  wire  position 
with  a  mean  pixel  value  at  near  by  positions.  It  is  difficult  to  remove  all  the  residual  trace  of 
the  wire  this  way  as  can  be  seen  from  the  following  image  display  especially  for  the  phase 
contrast  image  (Fig.  5(c)).  In  the  phase  contrast  image  (Fig.  5(c)),  the  small  residual  variation 
from  the  original  wire-track  really  got  enhanced.  The  linear  attenuation  coefficients  and  elec¬ 
tron  densities  for  the  50%  glandular-50%  adipose  breast  tissues  are  computed  from  the  tissue’s 
elemental  composition  and  the  interpolated  elemental  mass  attenuation  coefficients  from  the 
standard  reference  in  the  mammographic  radiation  dosimetry  [34,  35],  Moreover,  each  mass 
attenuation  coefficient  is  further  broken  down  to  two  components:  one  from  x-ray  photoelec¬ 
tric  absorption  and  coherent  scattering,  and  another  from  incoherent  scattering.  In  this  way, 
the  total  attenuation  is  partitioned  as  a  product  of  A“c  coh  and  as  defined  in  Eq.  (8)  above. 
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In  addition,  to  simulate  the  microcalcifications  in  breast,  four  small  ellipsoids  of  calcium  car¬ 
bonate  (CaC03)  are  embedded  in  the  specimen  model.  The  diameter  of  the  ellipsoids  can  be 
adjusted  in  simulating  different  size  of  the  microcalcifications.  In  the  following  simulations,  to 
test  phase-contrast  sensitivity,  we  set  the  diameters  of  the  four  ellipsoids  to  10,  5,  10,  and  5 
microns  in  x-ray  direction  and  300,  200,  300,  200  microns  in  detector  plane  respectively. 

The  attenuation  image  A 5  of  the  specimen  model  simulated  with  35.5  keV  x-ray,  and  its  two 
corresponding  partition  images  Aj^  and  Apecoh,  are  shown  in  Fig.  3.  For  a  comparison,  the 
profiles  of  A^  and  A^e  coh  along  the  line  passing  through  the  microcalcifications  are  shown  in 
Fig.  4  simulated  with  x-ray  energy  equals  18.5,  35.5  and  59.5  keV,  respectively.  We  can  see  that 
the  contribution  of  ApC  coh  to  the  total  attenuation  A5  gets  smaller  when  x-ray  energy  is  getting 
higher.  Especially,  when  x-ray  energy  equals  59.5  keV,  the  contribution  of  A^e  coh  for  the  soft 
tissue  can  almost  be  neglect,  as  is  expected. 

3.2.  Comparison  with  the  GS  Algorithm 


(a)  (6)  (c) 


Fig.  5.  Image  representation  of  the  inputs  generated  from  the  simulation  model  and  Fresnel 
propagation,  (a)  the  phase  map  <j>;  (b)  the  attenuation  map  Aq;  and  (c)  the  normalized  Fres¬ 
nel  propagated  phase  contrast  image  I  with  object  to  detector  distance  JO  =  26  in  (0.66  m). 

The  GS  algorithm  is  an  iterative  algorithm  for  phase  retrievals  from  a  pair  of  images  at 
two  planes  related  by  the  Fourier  transform.  For  details  readers  are  referred  to  [26],  The  GS 
algorithm  is  a  classical  algorithm  which  is  widely  used  in  electron  microscopy,  wave  front 
sensing,  astronomy,  crystallography,  and  many  other  fields  involving  phase  recovery  [27,  28, 
36], 

By  replacing  the  Fourier  Transform  in  the  GS  algorithm  with  the  Fresnel  transform  of 
Eq.  (11),  we  get  a  modified  version  of  the  GS  algorithm  extended  to  the  Fresnel  diffraction 
regime.  Our  previous  simulation  tests  showed  that  this  modified  GS  algorithm  converges  very 
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slow  for  object-detector  distance  Rs  ~  1  m,  and  converges  faster  for  larger  /G,  such  as  images 
acquired  at  synchrotron  beam  lines.  But  this  is  generally  not  applicable  to  the  held  of  clinical 
imaging,  due  to  the  physical  constraints  such  as  compact  sizes  of  hospital  rooms. 

In  our  simulation  tests,  we  compare  the  performance  of  APBA  and  that  of  the  GS  algorithm. 
The  photon  energy  of  the  point  x-ray  source  is  set  to  35.5  keV,  and  the  distances  from  the 
source  to  object  and  object  to  detector  are  set  to  R\  =  R2  =  26  inches  (0.66  m)  respectively. 
For  convenience,  the  incident  x-ray  intensity  Tm  at  R\  is  set  to  M2,  where  M  =  (R\  +  Rf)/R\ 
is  the  magnification  factor.  The  phase  map  (j)  and  the  attenuation  A5  are  generated  from  our 
phantom  model  for  35.5  keV  x-ray.  Fig.  5  shows  the  simulated  images  of  the  phase  map  (j).  the 
attenuation  image  A5  and  the  phase-contrast  image  I. 


(a)  (b)  (c) 

Fig.  6.  Comparison  of  the  performance  of  the  GS  algorithm  and  APBA.  (a)  plot  of  the 
accuracy  measures  with  respect  to  iteration  steps.  The  plot  with  solid  line  represents  the 
APBA.  The  one  with  dashed  line  represents  the  GS  algorithm;  (b)  recovered  phase  map 
with  the  GS  algorithm  after  100  iterations;  (c)  recovered  phase  map  with  APBA  after  100 
iterations. 

In  the  simulation  test,  the  iteration  for  our  attenuation-partition  based  algorithm  (APBA) 
and  the  GS  algorithm  are  performed  100  steps.  The  corresponding  recovered  phase  images  are 
shown  in  Fig.  6(b)  and  (c).  In  Fig.  6(a),  the  solid  line  represents  the  change  of  total  variation 
(TV)  of  the  retrieved  phase  using  attenuation-partition  based  algorithm  with  respect  to  the  it¬ 
eration  steps.  The  dashed  line  represents  the  change  of  TV  of  the  retrieved  phase  using  the 
(modified)  GS  algorithm  with  respect  to  iteration  steps.  We  can  see  that  the  converging  speed 
of  the  (modified)  GS  algorithm  is  much  slower  than  that  of  attenuation-partition  based  algo¬ 
rithm  (the  change  of  TV  of  APBA  from  step  1  to  step  100  is  1 .04 E  —  3,  almost  33  times  greater 
than  that,  3 . 172?  —  5,  of  the  GS  algorithm.)  In  addition,  from  the  visual  perception  point  of  view, 
we  see  that  the  phase  map  retrieved  with  the  attenuation-partition  based  algorithm  (APBA)  is 
much  better  than  that  retrieved  with  the  (modified)  GS  algorithm. 

The  main  difference  between  APBA  and  the  GS  algorithm  is  that  APBA  takes  the  advan¬ 
tage  of  the  phase-attenuation  correlation  property,  although  the  extent  of  this  correlation  is  not 
known  in  priori ,  but  the  GS  algorithm  does  not.  From  this  example  we  see  that  the  phase- 
attenuation  correlation  is  a  very  important  information  that  shouldn’t  be  neglected  in  the  algo¬ 
rithm  development  for  phase  retrieval. 

3.3.  Comparison  with  the  Transport  of  Intensity  (TIE)  algorithm 

We  have  mentioned  the  transport  of  intensity  equation  in  Eq.  (3)  in  Section  1.  Since  Teague 
proposed  the  TIE  algorithm  for  phase  retrieval  in  1983  [21],  numerous  phase  retrieval  algo¬ 
rithms  have  been  suggested  on  how  to  effectively  search  the  numerical  solution  of  the  TIE 
[4,  20,  25,  37,  38,  39,  40].  Among  the  methods  of  solving  the  above  TIE  for  the  phase  map, 
the  one  based  on  Fast  Fourier  Transform,  proposed  by  Nugent  et  al.  [4],  and  Allen  and  Oxley 
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in  [25],  is  the  most  widely  used.  In  the  form  of  pseudo-differential  operator,  the  solution  phase 
map  0  is  given  by 


<Hr) 


v[v-2(/-a2)] 


^5 


(14) 


for  the  given  normalized  phase-contrast  image  I  and  attenuation  image  A5.  Here  V  is  the  gradi¬ 
ent  operator,  and  V  2  is  the  inverse  Laplacian  operator. 

The  advantage  of  this  algorithm  is  that  it  does  not  require  the  boundary  information  in  solv¬ 
ing  the  TIE  (assuming  the  image  data  is  periodic);  it  is  a  deterministic  method  and  thus  the 
algorithm  is  fast  and  accurate  comparing  to  most  iterative  algorithms.  In  this  section  we  com¬ 
pare  the  performance  of  the  TIE  algorithm  and  that  of  APBA  for  two  kind  of  cases:  first  for  the 
ideal  case  without  any  noise  and  any  image  misalignment,  then  for  cases  simulating  practical 
situation  with  x-ray  imaging  noise  and  possible  image  misalignment.  In  these  simulation  tests, 
the  imaging  geometries  are  the  same  as  in  the  previous  subsection,  and  x-ray  energy  is  again 
35.5  keV  . 

For  the  ideal  case  without  any  noise  and  any  image  misalignment,  the  performance  compar¬ 
ison  results  are  shown  in  Fig.  7.  For  the  ideal  case  the  TIE  algorithm  is  accurate  both  in  TV 
measure  and  visual  perception.  APBA  can  also  achieves  this  accuracy  but  needs  1110  iteration 
steps  to  have  its  TV  measure  of  0.00215226,  an  error  smaller  than  that  of  0.00215269  with  the 
TIE  algorithm. 

However,  the  real  test  lies  in  the  performance  for  the  cases  simulating  practical  situation 
with  x-ray  imaging  noise  and  possible  image  misalignment.  Obviously  only  the  performance  in 
these  cases  really  matters  in  phase  contrast  imaging  applications,  especially  for  clinical  imaging 
applications  where  the  imposed  radiation  limits  dictates  existence  of  substantial  x-ray  quan¬ 
tum  noise  in  acquired  images.  Implemented  in  the  Fourier  space,  the  inverse  Laplacian  V  2 
in  Eq.  (14)  is  singular  at  zero  spatial  frequency.  This  singularity  will  amplify  the  noise  in  the 
images  and  result  in  failure  of  accurate  phase  retrieval  for  the  TIE  algorithm.  To  overcome  this 
difficulty,  some  kind  of  regularization  must  be  used.  The  most  widely  used  regularization  is 
Tikhonov  regularization,  which  replaces  V~2  by  |u|2/(|u|4  +  ff2),  for  some  “favorable”  con¬ 
stant  parameter  K,  called  the  Tikhonov  regularization  parameter.  In  this  regularization  scheme, 
the  singularity  is  regularized,  and  the  favorable  parameter  K  means  the  retrieved  phase  (j)K  corre¬ 
sponding  this  K  is  as  close  to  the  true  phase  0  as  possible.  Roughly  speaking,  the  regularization 


Fig.  7.  Comparison  of  APBA  and  the  TIE  algorithm  with  pure  data,  (a)  plot  of  the  accuracy 
measure  with  respect  to  iteration  steps.  The  plot  with  solid  line  represents  the  APBA.  The 
one  with  dashed  line  represents  the  TIE  algorithm.  It  needs  1110  steps  for  the  TV  measure, 
0.00215226,  of  APBA  to  achieve  to  the  TV  measure,  0.00215269,  of  TIE;  (b)  recovered 
phase  using  the  TIE  algorithm;  (c)  recovered  phase  using  APBA  after  1500  iteration  steps. 
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parameter  K  is  inversely  proportional  to  the  images  signal-noise  ratio.  Two  problems,  however, 
arise  to  this  regularization.  First,  the  true  phase  <j>  is  not  known  in  real  situations.  So  the  regu¬ 
larization  parameter  K  is  not  a  priori,  which  makes  it  difficult  in  practical  applications.  Second, 
this  Tikhonov  regularization  is  based  on  finding  a  stable  solution  to  Ax  =  y,  in  Hilbert  spaces 
X,  Y,  by  solving  the  minimization  problem 

xK  =  argmin  ||Ax  — y||y  +  k:||x|||  .  (15) 

xex 

It  is  Li  norm  dependent.  Since  the  proper  norm  for  image  data  is  the  total  variation  (TV)  norm 
[33],  a  favorable  solution  under  the  Tikhonov  regularization  principle  can  not  be  guaranteed 
a  best  solution  in  visual  perception.  Moreover,  for  relatively  noisy  acquired  images,  the  TIE- 
based  algorithm,  even  with  Tikhonov  regularization,  often  failed  to  retrieve  the  phase  maps[l]. 


(<*)  (e)  (/) 


Fig.  8.  Comparison  of  the  TIE  algorithm  and  APBA  with  noise  added,  (a)  True  phase  map 
</);  (b)  attenuation  map  A?;  (c)  the  normalized  phase  contrast  image  /;  (d)  recovered  phase 
map  with  the  TIE  algorithm,  no  Tikhonov  regularization  is  used;  (e)  recovered  phase  map 
with  the  TIE  algorithm  with  Tikhonov  regularization;  (f)  recovered  phase  map  with  APBA 
after  10  iteration  steps.  In  the  simulation,  the  acquired  data  is  assumed  to  have  a  level  of 
8i,  =  0.03%  detector  noise  and  one  pixel  misalignment  between  A^  and  /  horizontally. 

In  the  following,  we  will  compare  the  performance  of  APBA  and  the  TIE-algorithm  when  noise 
is  present.  In  the  practice  of  phase  retrieval,  there  are  generally  two  kinds  of  image  data  errors. 
One  is  the  noise  associated  with  image  acquisitions,  including  the  quantum  noise  of  x-ray  pho¬ 
ton  detections  and  detector  electron  noise.  We  assume  the  quantum  noise  dominates  as  is  the 
case  in  most  imaging  applications.  The  other  is  the  error  caused  by  the  misalignment  between 
the  attenuation  map  Aq  and  phase-contrast  image  I.  This  is  because  usually  the  attenuation 
image  and  corresponding  phase  contrast  image  are  generally  acquired  in  two  separate  x-ray 
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exposures.  The  misalignment  could  be  resulted  from  the  shift  or  tilting  of  the  object  or  detector 
between  the  exposures.  In  the  simulation,  we  associate  each  pixel  value  of  an  image  a  photon 
count  V,  so  that  P(i.j)  =c-N(i,j),  where  c  is  a  constant.  Assuming  the  noise  has  a  Poisson  dis¬ 
tribution,  with  variance  a2  =  N  at  each  pixel.  In  the  simulation  we  assign  a  background  noise 
level  for  each  simulated  image.  This  background  noise  level  is  defined  as  the  ratio  8b  =  (?b/Nb 
corresponding  to  the  direct  exposure  area  (where  A2  =  1)  outside  the  object  in  background.  The 
Poisson  statistics  dictates  that  Nb=  1  /§£  and  in  this  way  the  photon  count  N (i.  j)  can  be  deter¬ 
mined  accordingly  at  each  pixel.  Once  N(i,j)  is  determined,  the  statistical  errors  at  each  pixel 
can  be  assigned  using  a  computer  simulated  random  Poisson  distribution  generator  with  mean 
corresponding  to  the  photon  counts  N(i,j).  In  the  simulations  below,  the  background  noise 
level  8b  is  set  to  0.0003.  The  images  with  the  noise  added  are  shown  in  Fig.  8(b)  and  (c).  The 
quality  of  an  image  depends  not  only  on  the  noise  level  but  also  on  the  extent  of  image  contrast 
change.  With  the  assumption  P(i,j)  =  c -N(i,j)  above,  one  can  easily  see  that  8/,  =  (Jpb/Pb  is 
the  statistical  coefficient  of  variation  in  absorbed  photon  numbers  in  background.  8()i  =  OpJPoi, 
the  coefficient  of  variation  of  the  object  image,  on  the  other  hand,  is  the  structural  coefficient  of 
variation  of  the  sampled  image,  which  reflects  the  image’s  normalized  extent  of  image  contrast 
change.  We  will  use  the  ratio  8ol  /  8b  to  reflect  the  image  quality.  The  larger  the  ratio,  the  higher 
the  image  quality.  In  our  simulation  models,  when  8b  =  0.0003,  the  corresponding  ratio  8t)i/8b 
for  attenuation  A5  and  phase-contrast  image  /,  Fig.  8(b)  and  (c),  are  1.67  and  6.35  respectively. 
Because  of  the  phase-contrast  effect,  the  image  quality  of  the  phase-contrast  image  I  is  higher 
than  the  attenuation  image  A5  although  they  have  the  same  background  noise  level. 

Three  simulation  tests  are  performed  in  the  comparison:  easel:  assume  the  acquired  data  Aq 
and  I  has  noise  present  but  perfectly  aligned;  case  2:  assume  acquired  data  has  no  noise  but 
has  one  pixel  misalignment  horizontally;  case  3:  assume  acquired  data  has  combined  detector 
noise  as  well  as  one  pixel  misalignment.  One  bias  in  simulation  for  the  TIE-algorithm  should 
be  mentioned:  with  the  known  true  phase  value,  a  favorable  Tikhonov  regularization  parameter 
K  can  be  searched,  but  in  practice  this  search  is  hardly  feasible.  In  each  case  mentioned  above, 
three  phase  retrievals  are  performed:  1 .  using  the  TIE-algorithm  without  Tikhonov  regulariza¬ 
tion;  2.  using  the  TIE-algorithm  and  the  favorable  Tikhonov  regularization  parameter;  3.  using 
APBA  with  10  step  iteration.  The  Total  variation  (TV)  of  the  results  are  listed  in  Tab.  1  and 
the  recovered  phase  images  for  case  3  are  shown  in  Fig.  8(d)  -  (f).  The  results  using  Tikhonov 
regularization  are  better  than  those  not  using  Tikhonov  regularization  but  worse  than  those  us¬ 
ing  APBA.  The  influence  of  misalignment  to  APBA  is  little  but  disaster  to  the  TIE  algorithm. 
From  the  profile,  shown  in  Fig.  9,  along  a  line  passing  through  the  microcalcifications  we  can 
see  that  the  values  of  the  phase  recovered  from  APBA  is  very  close  to  the  true  phase  value,  but 
the  values  of  the  phase  recovered  from  the  TIE  algorithm  are  distorted. 

4.  Discussion  and  Conclusion 


Table  1 .  TV  comparison  of  the  TIE  algorithm  and  APBA.  In  the  table,  K  is  the  Tikhonov 
regularization  parameter,  A  represents  the  sampling  step-size  in  FT-space. 


TV  = 

Case  1 

Case  2 

Case  3 

TIE:  no  Tikhonov 

0.0406 

0.0808 

0.0886 

TIE:  with  Tikhonov 

K  =  0.7828A2 

0.0399 

K  =  4.4421 A2 

0.153 

K  =  4.1076A2 

0.0379 

APBA 

0.0283 

0.0086 

0.0290 
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Fig.  9.  Profiles,  along  a  line  passing  through  the  microcalcifications,  of  the  recovered  phase 
using  APBA,  the  solid  line,  and  using  the  TIE  algorithm  with  Tikhonov  regularization,  the 
dashed  line,  in  Case  3.  The  dash-dotted  line  is  the  true  phase. 


In  above  analysis,  we  assumed  that  the  x-ray  source  is  a  quasi-monochromatic  point  source.  In 
practice,  one  often  employs  conventional  incoherent  and  polychromatic  sources  such  as  x-ray 
tubes  for  imaging.  In  the  experiments  with  an  x-ray  tube  source,  the  previous  formula  Eqs.  (2- 
3)  of  the  phase  contrast  image  intensities  should  be  modified.  Since  in  the  APBA  method  the 
duality  transform  is  derived  based  on  Eq.  (3),  hence  the  Eq.  (6)  should  be  modified  accordingly. 
In  our  previous  works  we  have  studied  this  problem  [20].  With  the  Wigner  function  based  phase 
space  formalism,  we  have  proved  that  the  coherence  degree  of  a  finite-size  focal  spot  can  be 
accounted  for  by  the  optical  transfer  function  OTFq  jj  (h/M)  for  the  geometrical  unsharpness 
associated  with  the  finite-size  source  [20]: 
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where  /spot(4)  is  the  intensity  distribution  of  the  focal  spot.  We  found  the  generalized  TIE 
equation  with  an  x-ray  tube  source,  that  is,  the  x-ray  intensity  at  the  detector  is  given  by  [20]: 
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where  operator  (g>  denotes  the  convolution,  A 2  is  the  total  attenuation  of  the  imaged  object,  0  is 
the  spectrally  averaged  phase-shift  caused  by  the  object,  and  ( •  )  means  the  spectral  average. 
Compare  above  equation  with  Eq.  (3)  and  it  is  clear  that  the  TIE-based  phase  retrieval  method 
needs  only  two  modifications:  (i).  Fourier  space  de-convolution  of  the  measured  intensity  from 
OTFq  jj  (u/M),  (ii).  Replacing  wavelength  with  the  spectral-averaged  (A2)/(A).  In  the  same 
fashion,  the  duality  transform  D  defined  in  Eq.  (6)  should  be  modified  with  (i).  Fourier  space  de- 
convolution  of  the  measured  intensity  from  OTFqjj.{u/M)\  (ii).  A  replacement  of  the  parameter 
k  defined  in  Eq.(7)  with  the  spectral-average 


(ic) 


reRl 
2: zM 


#126430  -  $15.00  USD  Received  2  Apr  2010;  revised  1 1  Jun  2010;  accepted  9  Jul  2010;  published  15  Jul  2010 
(C)  2010  OS  A  19  July  2010/ Vol.  18,  No.  15  /  OPTICS  EXPRESS  16088 


Otherwise,  the  APBA  flow  chart  is  the  same  as  that  for  the  case  with  a  quasi-monochromatic 
point  source. 

Phase  retrieval  is  a  crucial  step  for  quantitative  imaging  such  as  reconstructing  the  3-D  dis¬ 
tribution  of  tissue  linear  attenuation  coefficients  and  refraction  indices.  However,  images  ac¬ 
quired  in  medical  applications  are  relatively  noisy,  due  to  the  radiation  dose  constraints,  with 
low  phase  contrast  effect,  due  to  physical  constraints  such  as  compact  sizes  of  hospital  rooms. 
An  phase  retrieval  algorithm  which  is  robust  to  noise  is  necessary  for  potential  medical  phase 
contrast  imaging.  In  [1],  the  authors  developed  an  algorithm,  called  attenuation-partition  based 
phase  retrieval  algorithm.  It  is  an  iterative  algorithm  which  takes  advantage  of  the  correlation 
between  the  attenuation  and  phase-shift.  The  phase  retrieval  results  from  experimental  images 
show  that  this  algorithm  is  fast  and  robust  [1],  In  this  work,  we  systematically  compared  the 
performance  of  this  algorithm  with  other  two  widely  used  phase  retrieval  algorithms,  namely 
the  Gerchberg-Saxton  (GS)  algorithm  and  the  Transport  of  Intensity  Equation  (TIE)  algorithm. 
The  systematical  comparison  is  conducted  by  analyzing  phase  retrieval  performances  with  a 
digital  breast  specimen  model.  We  show  that  the  proposed  algorithm  converges  faster  than  the 
GS  algorithm  in  the  Fresnel  diffraction  regime,  and  is  more  robust  against  image  noise  than 
the  TIE  algorithm.  These  results  suggest  the  significance  of  the  proposed  algorithm  for  future 
medical  applications  with  the  x-ray  phase  contrast  imaging  technique. 
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processed  using  two  algorithms  to  localize  events:  (i)  simple-threshold  and 
(ii)  weighted-centroid  methods.  The  processed  frames  were  added  to  give 
the  final  image.  Results:  (i)  Stationary  collimator:  The  image  showed  a 
grid-like  pattern  corresponding  to  the  septa  walls  of  the  collimator.  All  of 
the  hot  rods  in  the  phantom  can  be  identified,  (ii)  Moving  collimator:  The 
septal  pattern  is  blurred  out,  and  all  of  the  rods  can  be  clearly  identified. 
Conclusion:  The  same  MAF  can  be  used  in  both  nuclear-medicine  imaging 
and  x-ray  imaging  in  SPC  mode  for  dual-mode  imaging.  Currently  we  are 
limited  by  collimator  resolution  for  radionuclide  imaging.  Support  from: 
NIH  Grants  R01EB002873  and  R01EB008425. 
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Bedside  SPECT  Imaging  with  Pinhole  Collimation 

A  Cebula1  *,  D  Gilland1,  M  Studenski2,  (1)  University  of  Florida, 
Gainesville,  FL,  (2)  Thomas  Jefferson  University,  Philadelphia,  PA 

Purpose:  The  objective  of  this  study  is  to  evaluate  the  imaging  performance 
of  a  mobile  SPECT  system  with  a  pinhole  collimator  and  compare  the 
results  with  parallel  hole  collimators.  The  goal  is  to  obtain  both  planar  and 
tomographic  performance  measures  with  Tc-99m  and  F-18.  This  abstract 
presents  the  results  using  Tc-99m.  Method  and  Materials:  The  system 
utilizes  a  small  field  of  view  camera  with  pixilated  Nal  crystal  and  position- 
sensitive  photomultiplier  tubes.  The  pinhole  collimator  is  a  tungsten  knife 
edge  with  a  hole  diameter  of  3mm,  a  focal  length  of  12.5cm,  and  an 
acceptance  angle  of  90  degrees.  The  parallel  hole  collimators  have  been 
previously  described  along  with  their  imaging  performance  evaluation  [1]. 
The  following  performance  measures  with  the  pinhole  collimator  were 
obtained  and  are  presented  here:  count  rate  performance,  energy  resolution, 
flood  field  uniformity,  system  spatial  resolution,  system  sensitivity.  Results: 
The  maximum  count  rate  was  calculated  to  be  1.58x10s  cps  corresponding  to 
an  activity  of  161  pCi  at  30cm.  An  energy  spectrum  from  the  flood 
acquisition  demonstrated  an  energy  resolution  of  20%  FWHM. 
Magnification  corrected,  system  planar  spatial  resolution  was  1 .03  cm  with  a 
system  sensitivity  of  1.39  cps/pCi  at  a  source  to  collimator  distance  of  10 
cm.  Conclusion:  The  imaging  performance  of  a  mobile  SPECT  system 
with  a  pinhole  collimator  has  been  presented.  In  comparison  to  parallel 
hole  collimation,  the  pinhole  collimator  provides  superior  spatial  resolution 
and  sensitivity  at  distances  less  than  5cm.  However,  at  further  distances, 
pinhole  sensitivity  declines  while  parallel  hole  sensitivity  remains  relatively 
constant. 

Research  sponsored  by  the  United  States  Army  Medical  Research  and 
Material  Command  under  Award  No.  W81XWH-04- 1-0594 
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On  the  Development  of  On-Board  PET  with  Tomotherapy  Using  Open 
Dual  Ring  Geometry 

N  Darwish1  *,  T  Mackie2,  B  Thomadsen3,  C  Kao4,  (1)  University  of 
Wisconsin,  Madison,  WI,  (2)  University  of  Wisconsin  -  ADCL,  Madison, 
WI,  (3)  University  of  Wisconsin,  Madison,  WI,  (4)  Univ  Chicago,  Chicago, 
IL 

Purpose:  Positron  emission  tomography  (PET)  using  open  dual  ring 
geometry  was  investigated  as  an  on-board  system  for  functional  imaging  and 
PET  marker  tracking,  specifically  with  tomotherapy.  The  dual  ring  PET 
would  allow  measurement  of  both  inter  and  intra-ffactional  variation, 
facilitating  the  determination  of  treatment  uncertainties  and  improving  the 
delineation  of  tumor  volume  at  any  stage  in  the  radiation  treatment  delivery 
process.  This  study  demonstrates  the  field  of  view  (FOV)  of  various  PET 
ring  axial  gaps,  the  sensitivity  of  each  system  for  each  gap,  and  the  image 
quality  as  a  function  of  scan  time  and  activity  for  one  specified  design. 
Method  and  Materials:  Investigation  of  FOV  and  sensitivity  of  each  design 
was  accomplished  via  Monte  Carlo  simulations  with  GATE,  the  Geant4 
Application  for  Emission  Tomography.  Image  quality  was  investigated  with 
a  fully  3D  OSEM  iterative  reconstruction  algorithm  for  noiseless  data  to 
investigate  the  FOV  and  noisy  data  to  investigate  the  scan  time  at  a  specified 
activity.  Results:  Even  when  the  axial  gap  (G)  exceeded  the  axial  width  (W) 
for  each  open  ring  system  as  in  the  extreme  case  of  G  =  600  mm  (required  to 
avoid  actuators  holding  MLC’s)  the  axial  FOV  was  approximately  200  mm. 
A  continuous  axial  and  transaxial  FOV  (determined  by  PET  ring  diameter) 
of  360  mm  and  450  mm  respectively  resulted  from  G=W=  120mm. 
Reconstructed  NEMA  2001  phantom  images  for  activity  of  15  mCi/70kg 
and  scan  times  ranging  from  10  to  100  seconds  was  demonstrated  without 
the  effect  of  attenuation,  scatter,  and  random  events.  Conclusion:  For 


tomotherapy  treatment  200  mm  continuous  FOV  axially  was  sufficient  and 
was  achievable  with  all  gaps  simulated  in  the  range  from  120  mm  to  600 
mm.  Open  dual  ring  on-board  PET  has  the  potential  to  acquire  images 
within  the  time  limit  of  a  tomotherapy  treatment  session. 

MO-E-204C-04 

Improving  Robustness  of  Phase  Retrieval  in  X-Ray  Phase  Contrast 
Imaging 

XWu1  *,  AYan1,  H  Liu2,  (l)Univ  Alabama  Birmingham,  Birmingham, 
AL,  (2)  Univ  Oklahoma,  Norman,  OK 

Purpose:  The  phase  retrieval,  retrieving  tissue  phase  maps  from  the  x-ray 
phase  contrast  images,  is  an  important  task  in  x-ray  phase-sensitive  imaging. 
The  robustness  of  phase  retrieval  algorithms  is  of  critical  importance  for 
reducing  radiation  doses  in  clinical  applications.  We  show  that  the 
conventional  phase  retrieval  method  is  actually  unstable  against  the  quantum 
noise.  We  present  a  more  robust  phase  retrieval  method  that  we  developed. 
Method  and  Materials:  We  first  studied  the  phase  retrieval  by  means  of 
the  conventional  Transport  of  Intensity  Equation  (TIE)  method  for  phantom 
imaging  (such  as  a  nylon  air  bubble  wrap).  For  improvement  we  developed 
a  robust  phase  retrieval  method  based  on  our  notion  of  the  phase-attenuation 
duality  and  attenuation  partition,  and  derived  a  robust  iterative  phase 
retrieval  algorithm.  This  algorithm  had  been  applied  to  experimental  phase 
contrast  imaging  with  phantoms.  In  addition  the  phase  retrieval  accuracies 
had  been  analyzed  with  computer  simulations  of  phase  contrast  imaging  of  a 
digital  phantom  of  breast  tissue  characteristics.  Results:  We  found  that  the 
TIE-based  phase-retrieval  suffers  from  an  intrinsic  singularity.  The  phantom 
imaging  showed  that  the  TIE-based  phase  retrieval  method  is  unstable 
against  the  noise  and  misalignment  in  phase  contrast  imaging,  even  if  the 
Tikhonov  regularization  was  employed.  In  contrast,  our  method  based  on 
the  phase-attenuation  duality  and  attenuation  partition  is  singularity-free, 
and  is  robust  in  terms  of  the  phase-map  image  quality  and  phase-map 
accuracies  against  the  noise  and  misalignment,  as  is  verified  by  means  of 
experimental  phantom  imaging  and  computer  simulation.  Conclusion:  The 
conventional  TIE-based  phase  retrieval  method  is  unstable  against  noise  and 
image  misalignment.  Our  phase  retrieval  method  method  based  on  the 
phase-attenuation  duality  and  attenuation  partition  is  robust  against  the 
noise  and  image  misalignment. 

Research  supported  in  part  by  DoD  Breast  Cancer  Research  Program. 
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X-Ray  Luminescence  Computed  Tomography  Via  Selective  X-Ray 
Excitation 

G  Pratx*,  C  Carpenter,  C  Sun,  L  Xing,  Stanford  University  School  of 
Medicine,  Stanford,  CA 

Purpose:  X-ray  luminescence  computed  tomography  (XLCT)  is  proposed 
as  a  new  molecular  imaging  modality  for  imaging  X-ray-excitable 
phosphorescent  nanoparticles  three-dimensionally,  in  living  subjects.  Some 
of  these  nano-sized  particles  can  emit  near-infrared  (NIR)  light  when  excited 
with  X  rays  and  are  particularly  well  suited  for  in-vivo  biomedical  imaging 
because  the  signals  can  propagate  long  distances  in  tissue.  Method  and 
Materials:  The  imaging  mechanism  used  in  XLCT  consists  in  irradiating 
the  subject  using  a  sequence  of  programmed  X-ray  beams,  while  sensitive 
photo-detectors  measure  the  light  coming  out  of  the  subject.  By  restricting 
the  X-ray  excitation  to  a  single,  narrow  beam  of  radiation,  the  origin  of  the 
light  photons  can  be  inferred  regardless  of  where  these  photons  were 
detected,  and  how  many  times  they  scattered  in  tissue.  By  including  an  X- 
ray  detector  in  the  system,  anatomical  imaging  is  performed  simultaneously 
with  molecular  imaging  via  standard  X-ray  computed  tomography.  The 
molecular  and  anatomical  images  are  spatially  and  temporally  co-registered. 
Simulations  of  an  XLCT  system  were  performed  using  an  analytical  beam 
model.  A  preliminary  experiment  was  also  conducted  using  a  superficial 
treatment  beam  and  an  EM-CCD  camera.  Results:  Tracer  uptake  in  a  2  mm- 
diameter  target  can  be  detected  and  quantified  with  sub-picomolar 
sensitivity  using  less  than  1  cGy  of  radiation  dose,  a  result  that  makes  XLCT 
potentially  more  sensitive  than  PET,  currently  one  of  the  most  sensitive 
molecular  imaging  modalities.  Provided  sufficient  signal-to-noise  ratio,  the 
spatial  resolution  of  the  system  can  be  made  as  small  as  needed  by 
narrowing  the  beam  aperture.  In  particular,  1  mm  uniform  spatial  resolution 
was  achieved  for  a  1  mm-wide  X-ray  beam.  Images  reconstructed  from 
experimental  XLCT  measurements  showed  good  agreement  with  the 
simulation  model. 
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